#include "Eigen/Sparse"
#include "iostream"
#include <vector>
#include <cstdio>
#include "algorithm"
#include <Eigen/SparseCholesky>
 
using namespace std;
using namespace Eigen;
 
int main() {
    int i,j,m,n,l,a,d;
    int *I,*J,**J0,flag;
    double val, c;

    const int dim = 119425;        // 稀疏矩阵信息，这里偷懒假设事先知道！
    int nnz = 825737;

    /// 
    int maxcouple = 20;   // 每行最多有几个非零元？ 给一个适当大的数，请保证数组不越界的同时内存够用!
    char dummy;  // buffer for reading the data files

    FILE *fp, *fin;
    if ((fp = fopen("spmat.csv", "r"))==NULL)
    {
        puts("Fail!");
	exit(0);
    }

    SparseMatrix <double> A(dim, dim);
    for (int i = 0; i<nnz; i++)
    {
        fscanf(fp, "%d,%d,%lf", &a, &d, &c);
	A.insert(a-1,d-1) = c;
    }

    fclose(fp);

    A.makeCompressed();
    
    VectorXd rhs(dim);
    fin = fopen("rhs.csv", "r");
    if (fin != NULL)
	printf("右端项数据文件打开成功！！\n");
    else{
	printf("右端项数据文件打开失败，请检查！！\n");	return 0;
    }
    for(i = 0; i < dim; i++){
	fscanf(fin, "%le", &val);
	//	printf("rhs[%d] = %30.15f\n", i-1, val);
	rhs(i) = val;
    }
    fclose(fin);

    /**
     * step 4: Solving the linear system 
     */
    VectorXd solution;
    SimplicialLDLT<SparseMatrix<double> > solver;
    solver.compute(A);
    if(solver.info()!=Success) {
    // decomposition failed
    return -1;
    }
    solution = solver.solve(rhs);
    if(solver.info()!=Success) {
    // solving failed
    return -1;
    }

    cout<<"x=\n"<<solution<<"\n";    
    cout<<"error="<<(A*solution-rhs).norm()<<endl;
    return 0;
}
